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ABSTRACT 


A barotropic primitive equation model on a global tropi- 
cal grid using operational real data is used to test various 
boundary conditions and methods of initialization including 
numerical variational analysis. Comparisons of forecast 
accuracy are made between staggered and non-staggered grids. 

All prediction models produce better verification stat- 
istics than persistence, with the staggered grid verifying 
better than the full grid. This appears to be a result of 
two-gridlength noise being of greater magnitude in the full 
grid model. Less altering of individual synoptic systems 
occurs in the variational analysis initialization compared 
to two forms of the nonlinear balance equation. AS a re- 
sult, nondivergent variational analysis input fields to the 


prediction models result in the best verification values. 
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I. INTRODUCTION 


Numerical weather prediction in the tropics has devel- 
oped slower than prediction in the mid-latitudes for many 
reasons. The lack of data and less understanding of the 
actual physical processes governing the circulations in the 
tropics are probably the major deterrents to progress. In 
recent years, the increase of meteorological data and in- 
creased understanding of the dynamics in the tropics have 
Sparked interest in attempts at operational forecasting. 

Tropical weather prediction can be approached on either 
a global band grid or a global grid. Certain characteris- 
tics of tropical prediction make the global band advanta- 
geous over the global grid. In mid-latitudes, much of the 
baroclinic energy exchange occurs on the long wave scale 
while in the tropics, the input of energy on the cumulus 
scale is also an important contributor to the dynamics. 
This complicates the global prediction by making the tro- 
pics a subset of the entire globe. In mid-latitude regions 
of the Northern Hemisphere, data are routinely available at 
all the mandatory levels as well as significant levels. 

In the tropics, data availability is a different story. 
Upper level data near 300 mb are available from aircraft 
reports combined with scarce conventional rawinsondes. 
Gradient level and surface level data are also routinely 


available in numbers approaching those necessary for input 





to an operational model. However, the tropical regions are 
lacking in mid-tropospheric observations. The combination 
of a smaller scale of motion and the observational data 
available make data initialization in the tropics somewhat 
different than that in mid-latitudes. Reduction of the 
divergence equation to the balance equation, which is a 
generalized diagnostic equation used to relate wind and 
heights in mid-latitudes, is dependent upon a scale analysis 
in which the Rossby number is assumed less than one. This 
is not satisfied in the tropics however. Also, the optimum 
amount of divergence remaining in the initialized fields, 
which are used as input to a prediction model, may differ 
in the tropics from that in mid-latitudes. 

As a result of the above considerations, Fleet Numeri- 
cal Weather Central, Monterey now produces an operational 
wind and temperature analysis for a global band from 41S to 
60N at eight levels ranging from the surface to 100 mb. A 
more detailed description is presented later in this thesis. 
This analysis is a potential input to the primitive equation 
prediction model. Steinbruck (1971) conducted preliminary 
experiments on this global band grid using climatological 
data. He investigated the properties of various boundary 
conditions on the primitive equation prediction model and 
various diagnositic balancing schemes, including numerical 
variational analysis (NVA) developed by Dr. J. M. Lewis 
m972) for the global band. 


The purpose of this thesis is to extend this research 





to daily operational real data and to introduce a staggered 
grid system in an attempt to reduce memory and time require- 
ments on the computer. A deliberate attempt is made to 
avoid introducing averaging and diffusion in order to pro- 
duce as pure a forecast as possible and isolate true differ- 
ences between boundary conditions and the various techniques 


mor calculating initial data. 
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Bee SREMEELVESEQUATION PREDICTION MODEL 


A. BASIC MODEL 

A Peg Vayvex primitive equation model, developed by Har- 
rison (1970) and Elsberry and Harrison (1970, 1972), was 
used as a basis for the barotropic models used in these ex- 


periments. The equations for the ten-layer model are: 


ou -L(u) + fv - mee + Fx Ci) 
OV = -L(v) - fu - mo? + Fy (2) 
ot Oy 
23 = -L(6) + heat (3) 
a = -L(q) + moisture (4) 
1 000 
dt = “141 999) (5) 
Bm) 2. re) 
mee” fon (Gf oy Cg D1 (6) 
R R/C 

Ad = 0 C. A Caw P (7) 
where Wes). = 2 (2 us. | fe vs 0 

~ ox Gr?” oy Gr] + 35 (us) (8) 


and normal meteorological symbols are used with: 
A = finite difference 

= layer mean value 

1000 = 1000 mb surface 
m= map factor 


S = general scalar dependent variable 


im 





L = advective operator 

Linear stability is maintained by requiring the ratio 
of the space increment to the time increment be greater 
than the phase speed of the fastest waves, which are gravity 
waves. At 300 mb, and on the grid used, this time step is 
five minutes since the space increment near 60N is approxi- 
mately 75 nautical miles. 

All models use a forward time step in the first predic- 


1 
tion step: S = S° + At * function (S°). All subsequent 


prediction steps use a leapfrog scheme where: ghtl = gn-l 


+ 2 At * function (S®). 


B. NON-STAGGERED GRID BAROTROPIC PREDICTION MODEL 

Because the prediction models to be discussed here are 
barotropic, Eqns. (3), (4) and (7) are not needed and since 
all prediction is done at the 300 mb level, Eqn. (5) be- 


comes. 


99300 » 2 ‘ 
4 <h tee Gy ee Ey 


The frictional terms in Eqns. (1) and (2) are neglected. 

Input data to the models are the winds and geopotential 
at 300 mb. 

The finite difference formulas are based on Arakawa's 
(1966) flux form. This form of differencing attempts to 
preclude nonlinear computational instability by requiring 
that the advective terms conserve the mean square of any 


advected quantity. Therefore, mean square kinetic energy 
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is conserved by the nonlinear terms and nonlinear instabili- 
ty is avoided. All horizontal derivatives are centered two 
grid mesh derivatives. The horizontal advective terms for 
a general variable S with J the index in the y direction and 
Pethe index in the x direction are written: 

Su Su Su 


P) 
Ox i ) m m 
ieee J,I+1/2 J,I-1/2 


Ax 


Sy ,rtit Sg,1 





a 2 is - a a 
aca 
Ax 
die y oo eal 1 ou eM 
12 [ms My 7-1! 
Ax 


Note that this finite differencing only averages parameters 


over two grid points and one grid length. 


C. STAGGERED GRID BAROTROPIC PREDICTION MODEL 

Computation of space derivatives in finite difference 
form incorporates data from adjacent grid points. This 
prompted attempts at distributing the variables at alternat- 
ing grid points in an attempt to reduce computer time and 
memory. This approach was first attempted by Eliassen 
(1956) and subsequently has been used by many investigators 
in various forms. The prediction and conservation equa- 
tions on the staggered grid are the same as on the fea 
grid. The staggered grid distribution used in the baro- 


tropic prediction model involves u and v winds at alternat- 
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ing points, and w at diagonally alternating points as 
shown in Fig. 1. The Mintz-Arakawa general circulation 
model as described by Langlois and Kwok (1969) and an ocean 
circulation model by Haney (1971) also use this form of the 
staggered grid. This grid uses the least number of points 
possible and still retains the characteristic that the 
pressure gradient force is evaluated over the same space 
increment as in the non-staggered grid. Advective terms, 
however, involve a larger space increment. 

The horizontal advective terms for a wind variable X 


are written: 


¢_ Xu, = (Xu) _ Xu 
ee) = (—) Ge) 
x=J,I x=J,It1/2 x=J,1I-1/2 


2Ax 


Xy,itit Xy 7 1 oy +3 
ea =~) 1- 
2 e+1 jee 
2Ax 

AES ati i re-35 ru } 

2 2 m 3,1 fe 1 

2Ax 
Refer to Fig. 1 for indexing. In a similar manner, the 


horizontal advection of geopotential (%) takes the form: 





3 X X X 
e- (24) 2 (&) - 
O=J,I O=j,1+1/2 0=J,1-1/2 
2Ax 
‘ 3,14 i °IeI, res ares 
2 Z Jee Tm J~1,I1 
Zi 





STAGGERED GRID 


J+1 


Be rk 





FORECASTING GRIDS 
GLOBAL BAND: 49 x 146 (COLUMN 145 = 1, COLUMN 146 = 2) 
STAGGERED GRID: U,V = 25 x 74 (COLUMN 73 = 1, COLUMN 74 = 2) 
#,0 = 26 x 74 (COLUMN 73 = 1, COLUMN 74 = 2) 
STAGGERED DATA EXTRACTION 
U,V (321,25, 122,73) = U,V (2*J-1, 2*I-1) 
4,0 (J=2,25, 122,73) = €,w (2*J-2, 2*I-2) 
ROW 1 AND 26 DETERMINED BY BOUNDARY CONDITIONS. 


Fig. 1. Staggered grid description 
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aa os iba 1 
ee ee aa x X 
[ J [ Se sees =) 


J-1,1-1) 
2Ax 

The vertical motion used in vertical advection of wind 

speed is obtained by averaging the four nearest diagonal 

values of Ww. 

As in the full grid, this form of finite differencing 
averages parameters over only two points but results in 
differencing over two full grid lengths. 

Vertical motion is computed from the continuity equa- 
tion using winds that are obtained by averaging predicted 


winds over two grid increments located diagonally from the 


vertical motion gridpoint. 


D. BOUNDARY CONDITIONS 

The computationally stable integration of the primitive 
equations requires a set of boundary conditions that are 
consistent with the finite difference scheme. 

Flux through the top of the domain is eliminated by set- 
Eang &© = O at 100 mb. 

Two sets of horizontal boundary conditions were tested 
in the prediction model: 

1. No-flux boundaries 

Elsberry and Harrison (1971) discuss the boundary 

conditions that were applicable to the basic ten-level model 
and which are used in one version of the barotropic model. 
In the full grid, no flux of mass or energy is allowed at 


the northern boundary by setting: 


16 








V = _ U = J = 
ee gees 


where J+tl is the northern-most row in the grid. The south- 
ern boundary is treated in a similar manner. These boundary 
conditions result in the model boundary being a fictitious 
grid row between the two outer rows of input data across 
which there is no flux of mass or energy. East~-west bound- 
ary conditions are not needed since the global band is con- 
tinuous. 

The u boundary conditions on the staggered grid are 
Similar to those on the full grid. In order to conserve 


kinetic energy on the staggered grid ® 54. = 5 where J+l is 


1 
the northern-most row of the staggered geopotential values. 
To satisfy the no~flux criterion, the outer row of v com- 
ponent wind must be identically zero. This specification 
of boundary conditions results in the model boundary being 
coincident with the outer-most row of staggered grid wind 
components. 
2. Restoration boundaries 

Kesel and Winninghoff (1970) describe a constant 
flux, restoration boundary technique. This technique is 
used with initial data from the nonlinear balance equation 
which has consistent boundary conditions, and balancing 
produced by a numerical variational analysis scheme. Both 
of these techniques are later described in detail. The 


restoration technique involves restoring the newly predicted 


values near the boundaries toward their value at the previous 


IEE 





time step with a specified restoration coefficient. In the 
full grid prediction, these coefficients are a linear varia- 
momeror 1,0, 0:8, 0.6, 0.4, 0.2, 0.0. This results in a 
persistence forecast on the outer boundary, pure dynamical 
forecast six rows in and a linear combination between. Note 
that the data boundary is the model boundary. [In the 
Staggered grid a similar approach is used with coefficients 
ome .0, 0.6667, 0.333, 0.0. In the staggered grid, a 
fictitious geopotential value is needed one grid length 
outside the boundary of the available data so a linear 
extrapolation on the input data geopotential gradient is 


used: 
ieee ais = Ne 


This results in the model geopotential boundary being one 
row outside the input data boundary and the model wind 


boundary coinciding with the data boundary. 
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ITIL. DATA DESCRIPTION 


First test experiments with the prognostic and diagnos- 
tic models used a climatological input produced by Jenne 
(1969) of the National Center for Atmospheric Research. In 
particular, these experiments used 300 mb January climato- 
logy for the area of the global band. Although this data 
lacks the small-scale features which are important in day- 
to-day forecasting, the strong interaction between the tro- 
pics and the middle latitudes through the subtropical jet 
is present. In addition, boundary condition effects and 
development of near-equatorial circulations can be tested. 
The 300 mb level was chosen so that comparable later experi- 
ments using real data would allow a maximum of actual 
observations. 

The real data utilized in later experiments are from 
the Fleet Numerical Weather Central global band analysis 
developed by Grayson (1971). This analysis covers the 
entire global band from 60N to 41S. The grid dimensions 
are 49 X 144 on a Mercator secant projection map true at 
fgeo° latitude. The mesh length is 2.5° longitude (=150 
mimi.) at the equator and reducing to 75 n.mi. at 60N. This 
grid size has the capability of depicting relatively small 
scale synoptic events, provided the observations are suffi- 
ciently dense. Grayson (1971) used a modified version of 


the successive approximations technique to analyze sea 


iy 





level pressure and surface wind. This approach has subse- 
quently been extended to the analysis of seven levels of 
upper level winds, but it should be noted that no vertical 
consistency check is applied. In an attempt to accommodate 
the large variation of data coverage and the various scales 
of motion in the wind analysis, the area influenced by each 
observation, in addition to being dependent on the scan 
number, is dependent upon the reported wind direction (if 
greater than 13 knots). This approach was suggested by 
Endlich and Mancuso (1968). 

The previous six hour old analysis is used as a first 
guess field. 450 reports per analysis is a typical number 
of reports at 300 mb with approximately the following lati- 


tudinal distribution: 


45N - 60N = 372 
30N - 45N = 33% 
so Gee 1 0) 04 
158 - 15N = 62% 
S10 See Se — 
41S - 308 = 2% 


No upper air geopotential fields are analysed by the 
global band analysis program. Therefore, a first guess 
field is produced using a combination of the current FNWC 
hemispheric analysis north of 27N and climatology south of 
17N with a linear combination in between. 

Since the real winds and climatological heights along 


the southern boundary of the global grid are not geostrophi- 
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cally consistent, the height values are altered along the 
boundary. Geostrophic gradients appropriate to the 
analysed v component of the wind are added to the average 
climatological height along the southern wall. This 
approach is applied to all diagnostic models except no-flux 


boundary conditions where it is not applicable. 
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Woe litt LAL LZALLON METHODS 

A. NONLINEAR BALANCE EQUATION WITH NO-FLUX BOUNDARY 

CONDITIONS 

Elsberry and Harrison (1971) utilized the nonlinear 
balance equation and applied no-flux boundary conditions to 
formulate a diagnostic model which matches the no-flux ver- 
Sion of the prediction model. The nondivergent part of the 
wind is obtained from ¢ = vy where € is the relative vorti- 
city obtained from the observed u and v wind components. 
No-flux boundary conditions imply that the north and south 
boundaries are streamlines with no v component wind across 
them. Conservation of kinetic energy requires that the 
two outer rows of ¢ at the boundaries be set equal. [In 
order to avoid setting up gravity waves propagating along 
the boundary, the geopotential along the entire boundary 


is set equal to the average value. y 


v 


rss is set equal to 


Mee th -uAy where u is the mean zonal wind component 
averaged over the grid and Ay is the distance from the north 
to south boundary. The natural cyclic east-west boundary 
conditions alleviate the need for artificial boundaries in 
this direction. 

The nondivergent wind is calculated from the relaxed jp 
field using v = k x Vw “eo -Ttluxscomditions are foreed py 
making v(49,I) = v(48,1I) m(49)/m(48) and similarly at the 


southern wall where m is the map factor. Then, the geo- 
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2 
potential field is generated by relaxing V $ using the non- 


linear balance equation: 


2 


2 
= Ve 97 oy oy 
Vo VE Ve oy ote 2 (ax ay) 


One should note that large deviations in the geopotential 
from the input geopotential are to be expected due to the 
zonal averaging of the boundary values. 
B. NONLINEAR BALANCE EQUATION WITH RESTORATION BOUNDARY 

CONDITIONS 

Kesel and Winninghoff (1970) describe a restoration 
technique which eliminates the need for altering the bound- 
ary geopotential field and appears to control false reflec- 
tion of physical and computational modes at the boundary. 
Wind and height values of the input data are not altered 
on the north-south boundaries in contrast to the no-flux 
method. The technique involves restoring after each time 
step the values of the parameters near the walls toward 
the value in the previous time step with a specified restor- 
ation coefficient. This results in a purely dynamic fore- 
cast in the interior, a persistence forecast on the boundary 
and a combination in-between. The technique results in the 
boundaries acting as an energy sponge for externally pro- 
pagating meteorological and gravity waves. 

When working with climatological fields for input data, 
the southern boundary of W was determined similar to the no~- 
flux case by using a column average of the zonal wind conm- 


ponent. This technique applied to real data, however, pro- 
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duced unrealistic meridional velocity components along the 
southern wall due to the large and rapid variations in the 
east-west direction of the average zonal wind component. 

To correct this problem, the W's calculated as above are 

averaged. The southern ~ value is set to this value which 
results in a north-south WU gradient consistent with the 

mean zonal wind. Upon this average value at the southern 
boundary, variations are superimposed to reflect east-west 
gradients which are consistent with the analysed v compon- 
ents on the southern wall. The remainder of the balancing 


is identical to the no-flux approach. 


C. NUMERICAL VARIATIONAL ANALYSIS 
Lewis (1972) has adapted Sasaki's (1958) variational 
method to obtain dynamically consistent initial fields from 
input objective operational analyses. The fields are ad- 
justed so they satisfy the equations of motion and minimize 
the functional: 
2 2 2 


= ay oe pees = ye PUN a 
I J i fa¢u-u) + A(v-v) + B(Z-Z) + a (—) + a(x) Jas (9) 


meen tu, v, 2 being the input analysis. 4d, B are weighting 

functions specified beforehand which can be varied in res- 

ponse to the density and reliability of the input observa- 

tions and which determine how closely the desired fields u, 
v, and z fit the input analysis. The dynamic weight a is a 
Similar function which controls the horizontal acceleration 
(calculated from the momentum equations given below) and 


directly couples the wind and height information via the 
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assumed dynamical relations. 


Pieeee Oo Che Cie ce 
fe. wee Vv Vier ay g Vv Vv 258 


Q 

rr 

Qe 
Md 


The use of these equations as dynamical constraints require 
the simultaneous solution of the three equations for u, v, 
and z (so-called Euler equations). To simplify the Euler 
equations, the inertial terms in the dynamical constraints 


are modified: 


ou, ~ p22 _~ Ys yu + fv Ol ea pe 2 U5 SS Ro 





As in the restoration balancing, the climatological southern 
height values are averaged and altered to make them geo- 
strophically consistent with the v component winds before 


the NVA analysis is made on the fields. 
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VY. DIAGNOSTIC RESULTS 


Diagnostic experiments were aimed at three objectives: 
compare NVA results with various forms of the nonlinear 
balance equation; determine NVA weighting constants which 
most advantageously initialize the prediction model; and 
test various methods of limiting the divergence in the NVA 
output. 

The first objective was to compare the results from the 
variational analysis approach against those using the non- 
linear balance equation with no~-flux and restoration bound- 
ary conditions. One should first note the basic difference 
in solution approach between the two methods. The non- 
linear balance alters the geopotential field based on the 
nondivergent wind input from the streamfunction solution, 
while the variational analysis simultaneously alters both 
the winds and geopotential. 

Area average wind changes listed in the tables, and 
wind changes analysed in the charts to follow are obtained 
by subtracting the input field from the balanced field. 
Therefore, a positive value means an increase of westerly 
winds in the u component and an increase in southerly winds 
in the v component. Area average modulus change is the 
average of the absolute value of the individual wind 
changes which gives a measure of the magnitude of the 


average change at each grid point. The values given are an 
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average of the u and vy component changes which were found 
to have similar values. The RMS change is the square root 
of the average of the squares of the balancing changes. 
RMS u and v components are averaged since similar values 
were obtained. 

Fig. 2 shows a typical v wind component input to the 
objective analysis programs. Fig. 3 shows the v component 
changes produced by the nonlinear balance solution with 
restoration boundary conditions. The nonlinear balance 
approach tends to reduce the absolute magnitude of the wind 
maxima in both the u and v component wind with the centers 
of maximum change being of nearly equal magnitude. This 
wind change represents the divergent component of the wind. 
As shown in Table 1, the modulus and root-mean-square change 
in the no-flux approach are larger than in the restoration 
approach due to the crude boundary approximations. The 
non-realistic boundary conditions in the no-flux approach 
‘only produce significant differences from the restoration 
boundaries within about five rows of the boundary. The wind 
speed changes produced by the variational analysis approach, 
as displayed in Fig. 4, show smaller magnitudes in a larger 
scale pattern. Generally, the tendency is to reduce the 
Wind maxima and increase the intervening wind minima. Table 
1 shows an average modulus wind change of 0.60 m/sec in the 
pure NVA analysis as compared to a modulus change approxi- 


mately four times as large for the nonlinear balance method. 


Aay | 
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Fig. 2. 300 mb Global band v component wind analysis from 


60E to 150W for 08/1200Z January 1972. Units = 
m/sec. 





Fig. 3. V component wind speed changes resulting from non- 


linear balance solution with restoration boundary 


conditions. Positive (+) value indicates increase 
in southerly wind speed. Units = m/sec. 
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TABLE 1 


Balancing Area Area RMS 
Type Average Average Wind 
Wind Modulus Change 
Change Wind (m/sec) 


(m/sec) Change 


(m/sec) 


We Resulting Field Nondivergent 


NVA 


Wie 7,25) u=.50 
Made v=.2/7 


Nondivergent 


No-flux u=- 


Balance v= 
Equation 


Restoration 


Balance u=- 


Equation v= 


26 
223 


03 
~24 


prope “ised ocd 
84 4.12 
89 Zaye 


II. Resulting Field Divergent 


NVA u= 
mle7,25) v= 


go 


2025 


60 iO 
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METHODS OF BALANCING 


Area 
Average 
Height 
Change 
(m) 


-22.0 


106.7 


-22.0 


Area 


Average 
Modulus 
Height 
Change 
(m) 


ules se 


Se 
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0 


RMS 
Height 
Change 
(m) 


Sal) 


io 


174.4 


Sy 38 











Fig. 4. V component wind speed changes resulting from NVA 
inGenralezatton propram (1,7,25). Positive (+) 
value indicates increase in southerly wind speed. 
Units = m/sec. 





4 Bae — 





— 


—_ = == we 
i ee ee ee ee —™~2 aw. 








_— oe ean { bo) 
t— — : 4 ‘ ; oe 
Ff cS EE 
Gui = 

— eee me Pe eal O08 /izz 
yp 78 ge 
-——_ — my —_— oon oe —eae oo 

i I et ey ee 

= Ce ee Ra 





Fig. 5. 300 mb Global band height field from 60E to 150W 
tous, P2002 Janwanyerg/2. Units = decameters. 
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Fig. 5 shows the 08/1200Z January 1972 global band 
height field used in these experiments. In evaluating 
height changes produced by the various methods, ae must 
remember that the geopotential field was obtained from the 
derived nondivergent wind analysis. In addition, approxi- 
mately two-thirds of the height field, that south of 1/7N, 
is a climatology field. Keeping this in mind, one can 
see from Table 1 that the NVA analysis alters the heights 
only about one-quarter as much as the nonlinear balance 
approach. In the region with real height data, the NVA 
analysis tends to fill troughs and reduce ridge heights in 
order to achieve balance. These changes appear to be dis- 
placed toward the downwind side of the synoptic feature in 
the mid-latitudes. The scarcity of data in the southern 
hemisphere tends to produce isolated wind maxima and the 
NVA introduces height gradients to partially balance these 
winds. Fig. 6 shows that the NVA height changes necessary 
to generate gradients consistent with the winds are both 
positive and negative, and the changes are similar in 
scale to the wind maxima themselves. By contrast, the 
changes in the two cases using the nonlinear balance equa- 
tion approach are on a much larger scale as shown for the 
no-flux case in Fig. 7. The patterns are on the order of 
the long waves that are present, and the changes tend to 
deepen the troughs by building a ridge downwind of the 
trough axis. Unlike the NVA results, all height changes are 


positive in the northern hemisphere and are about four times 


Sal 








Fig. 6. Height changes resulting from NVA initialization 
peoeramn (ie 7.25). Units = m. 











Fig. 7. Height changes resulting from nonlinear balance 
solution with no-flux boundary conditions. 
Units = m. 
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as large as those in the NVA. Area average changes are 
slightly larger in the no-flux boundary condition than in 
the restoration case, since the boundary heights are con- 
stant. 

Early development experiments made use of climatologi- 
cal data as mentioned before. One would expect these fields 
to be composed of long wave features with none of the short 
wave systems apparent in daily analyses. This is evident 
in the wind and height changes necessary to balance the 
real data fields. The wind changes due to initialization 
of the climatological data are one-fourth to one-half of 
those for real data. 

The second objective of the diagnostic experiments was 
to develop weighting values for the variational analysis 
scheme which produced fields with sufficient balancing to 
successfully initialize the prognostic program. Table 2 
shows the results of this experimentation. Since only a 
wind analysis is available in the tropics, the NVA analysis 
should be forced toward the wind values with little weight 
on the height values. The dynamical constraint should be 
large enough to assure compatibility as an initial state 
for the prediction model. One must pay the price of addi- 
tional computer time if the analysis is to be forced toward 
exact balance, or small deviations from the input wind 
fields, as shown in Fig. 8. An arbitrary root-mean-square 
deviation in the wind field of 1.0 m/sec was set as a 


desirable limit on changes allowed from the input value. 
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TABLE 2 


VARIATIONAL ANALYSIS WEIGHTS 


a@ = winds 
B = geopotential 
a = dynamics 
Balancing Area Area RMS Area Area RMS 
Weights Average Average Wind Average Average Height 
B, a, a Wind Modulus Change Height Modulus Change 
Change Wind (m/sec) Change Height (m) 
(m/sec) Change (m) Change 
(m/sec) (m) 
i 
ies 25 u=1.43 1.40 2 -16.5 29.0 3D. 6 
v= .13 
Le u= .52 ~63 0 -22.1 347..0 37.0 
v= .025 
ie 6,25 u= .26 ee) .60 -25.2 34.8 Za 
v= .03 
2,209 u= .17 meee 741 -26.5 351.16 437.7 
v= .02 
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5 10 15 20 25 
RMS WIND CHANGE (m/sec) 


NVA wind weighting factor (solid), and number of 
iterations required for convergence (dashed) for 
a RMS wind difference between final and initial 
fields. The values of pressure and dynamical 
weighting factors are held fixed. 
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This resulted in an G@ value of approximately 7.0 when con- 
bined with B = 1.0 and a= 25.0 where 8 ima Seeneao ht 
changes, a specifies the degree of horizontal acceleration 
in the wind field and a@ limits the wind changes as speci- 
fied in Eqn. 9. One should note that limiting the wind 
deviations to very small changes causes a rapid increase in 
the number of iterations required for convergence. Speci- 
fying the coefficients of 1,7,25 required approximately 23 
seconds of computer time for each iteration of the uncon- 
verged points. This results in a total run time of about 
32 minutes which is extremely long for an operational analy- 
Sis. In contrast, the nonlinear balance approach requires 
approximately eight minutes of computer time. 

It was discovered in the experiments with climatologi- 
cal data that the wind fields resulting from the variational 
analysis program produced excessive two-gridlength gravity 
waves and excessive vertical motions. Eventually, this 
noise led to computational instability. The NVA procedure 
Satisfies the dynamical constraints on the fields by limit- 
ing the magnitude of horizontal accelerations, which are ob- 
tained from the wind tendency equations. There is no direct 
control on the amount of divergence allowed. Therefore, the 
third objective of the initialization experiments was to 
evaluate various methods of limiting divergence in the NVA 
Output field. Two alternatives existed; make the input 
wind fields to the NVA analysis nondivergent and produce a 


divergent output; or make the output of the NVA program non- 
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divergent by using only the nondivergent component. However 
this postadjustment of the wind will create some inconsist- 
ency ple. the height field that was produced by the NVA 
scheme. 

A FNWC subroutine, which uses version III of the method 
described by Hawkins and Rosenthal (1965), was used to conm- 
pute the nondivergent wind. It should be noted that this 
approach does not retain the original u and v winds on the 
boundary as boundary conditions. Nevertheless, changes are 
minimal on the boundaries. 

Table 3 shows the results of these experiments. The 
patterns of wind changes necessary to make the fields non- 
divergent show patterns similar to those of the nonlinear 
balance equation, however, on a smaller scale, as shown in 
Fig. 9. The individual wind maxima are altered similar to 
the scale af the NVA changes. One can see that the NVA pro- 
gram changes the winds approximately the same amount whether 
the input is divergent or nondivergent; however, the height 
changes necessary with the nondivergent input are larger. 
The modulus and root-mean-square wind changes necessary to 
make the fields nondivergent are about two and one-half 
times larger than changes in the NVA analysis. Neverthe- 
less, the total change involved in both methods are very 
nearly equal. It should be pointed out that a property of 
a nondivergent field is that a cyclic row average of the v 
component wind around the global band equals zero. This 


implies no net mass transport in the meridional direction. 
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TABLE 3 


VARIATIONAL ANALYSIS METHODS 


Method Area Area 
Average Average 
Wind Modulus 
Change Wind 
(m/sec) Change 
(m/sec) 
Method 1 
ae 
Divergent u=.25 5 pe 
enpiu t v=.2/ 
Made 
Nondivergent 
b. 
Nondivergent 
Input To u=.50 52) 
NVA v=-.02 
or.7,25) 
Total u=./5 ie 0 
Change v=.25 
kok Kk kK kK KK Kk OK KK 
Method 2 
Ae 
Divergent 
Input To u=.52 -60 
NVA v=.025 
57,2.) 
b. 
NVA 
mie? ,25)° u=—-.02 1.67 
Output v=.25 
Made 
Nondivergent 
Total u=.50 Jb teal 
Change v=.2/ 


RMS Area 
Wind Average 
Change Height 
(m/sec) Change 
(m) 

Zod None 

iL 4 (VG. -31.4 
Zoe -~31.4 
kok kK kK kK Kk KK 
1.10 -22.0 
2344 None 
226) -~22.0 
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Fig. 9. V component wind speed changes resulting from 
nondivergent subroutine. Positive (+) value 
indicates increase in southerly wind speed. 
Units = m/sec. 
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Vilweet Ke OE CLTEONS RESULTS 


A. BOUNDARY CONDITIONS 
Boundary conditions on a global band pose problems dif- 
ferent than those in a hemispheric model. The global band 
boundaries are in the mid-latitudes where strong meridional 
flow causes strong cross-boundary fluxes. A hemispheric 
model avoids these difficulties, since boundaries are loca- 
ted in the tropics where pressure gradients are relatively 
small and large scale waves are not commonly present. The 
major problem in the prediction models caused by the large 
cross-boundary flow is the generation of gravity waves at 
and near the boundaries. These waves contaminate the pre- 
diction of the meteorological waves and may lead to computa- 
tional instability. Prediction on the full grid with no- 
flux boundary conditions produce only small amplitude gravi- 
ty waves. Restoration boundary conditions with compatible 
nonlinear balance equation input to the prediction model 
produce strong two-gridlength gravity waves, which require 
filtering to obtain a stable 24-hour forecast. This filter- 
ing is accomplished using a scheme devised by Shapiro (1971) 
in which the variable's value at J,I is altered by 1/8 
S 5 4 


S ) where S is a 


CS 54417 joe ee ae eee a 2 


general variable. In order to obtain a 24-hour forecast, it 
waS necessary to apply the filter once every four hours. 


Figs. 10 and 11 show the gravity wave effect in the RMS 6 
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Grid Row for full grid restoration 


Fig. ll. Wrpms vs. 
Units = mb/sec. 


prediction. 
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and RMS w values respectively. The 2-d nature of the gravi- 
ty waves is evident in both fields. The cross-boundary 

flow and resulting gravity waves produce much larger RMS 6 
and Ww values near the northern wall which increases in 
amplitude with time. NVA input to the full prediction model 
resulted in the most severe 2-d waves, as might be expected 
since the dynamic constraint is not as strong as in the 
nonlinear balance approach. It was necessary to apply the 
filter twice at four-hour intervals to produce a 24-hour 
forecast. The staggered grid predictions did not require 
any filtering to produce a stable 24-hour forecast. RMS 9 
and RMS w for the staggered grid restoration prediction, 
shown in Figs. 12 and 13, are much better behaved than in 
the full grid prediction. At t=0, the w curve in Fig. 13 
arises from extracting alternate u and v component winds 
from the balanced solution obtained on the full grid. There- 
fore, the staggered field contains divergence while the full 
grid field is nondivergent and contains no vertical motion 
at t=0. The vertical motions in the staggered prediction 
Maintain the same size W as time increases. The amplitude 
of the two-gridlength waves is not as great and values at 
the northern boundary do not increase in amplitude as rapid- 
ly with time as in the full grid prediction models. This 

is to be expected since the finite differencing involves cal- 
culations over two gridlengths of the full grid and results 
in the use of points which have a similar phase angle with 


respect to the 2-d waves. 
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Fig. 12. dpms vs. Grid Row for staggered grid restoration 
prediction. Units = nm. 
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B. PREDICTION VERIFICATION METHODS 
The output fields were smoothed with a filter which 
altered the parameter value at J,I by 1/8 (S547 ,178 5-101 


ge I+1tS -4 § )}. Verification of the prognosis 


Wil J,I 
was based on these filtered fields. Also, to avoid un- 
realistic fields produced by the restoration process near 
the boundaries, and fictitious boundary specification in 
the objective analysis, verification statistics were only 
computed in the interior region ners a true dynamical 
forecast was computed. This excluded the outer six rows 
in the full grid prediction and the outer four rows in the 
staggered grid. Verification criteria for predictive 
models are numerous and by no means obvious. The area 
average modulus wind component differences and the RMS dif- 
ferences were chosen since accuracy in the wind analysis is 
more desirable in the tropics than are good height values. 
A single numerical value which was used to combine the 


various measures of forecast ability was defined as in the 


following formula: 


Forecast Value = (| CUR [lee sare. 
(|v | ; -|v| ) + (u - u ) 
t e 
pecraer © Verif RMS RMS 
Persist Verif 
ae (Vv Ny ) (10) 
RMS RMS 
Persist Verif 


where |u| area average modulus error 
u = n 
RMS root mean square error 


Persistence forecast error is defined as the comparable 
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balanced field at oe hours minus the comparable balanced 
field at toe Lavoie and Wiederanders (1968) have pointed 
out that in the tropics a persistence forecast plus a varia- 
ble amount of climatology is still difficult to beat with 
any objective prediction technique. 

Typical magnitude persistence forecast error values are 
shown in Table 4. 

Verifying forecast error is equal to the predicted 
values at t +24 hours minus the verifying balanced field at 
t +24 hours. Therefore, if the forecast value is positive, 
the verification is better than persistence with increasing 
magnitudes denoting better verification skill. 

A negative value implies less skill than a persistence 
forecast. Less emphasis was placed on height verification 
since these fields are derived rather than analysed. How- 
ever, one measure of prediction accuracy is the ratio of 
RMS persistence height change to the RMS verification 
height error. Numbers less than one imply a persistence 
forecast is more accurate. Values of this ratio greater 
than one imply forecast skill compared to persistence. 
These values are shown in Tables 5 and 6 in parentheses. 

The number of passes through the predicted fields with 
the smoother had a definite effect on the verification 
statistics. A typical example with the forecasts with NVA 
input is shown in Table 4. The modulus and RMS verification 
values for the v component are better than persistence, 


regardless of the number of filter passes. Successive 
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TABLE 4 


VERDPECATION STATISTICS As FUNCTION OF FILTER PASSES 


Forecast: Non-divergent NVA input (1,7,25) 


Number Area Area RMS RMS Area RMS 
of Average Average u Vv Average Height 

Passes Modulus Modulus Component Component Modulus Differ- 
u Vv Wind Wind Height ence 
Compon- Compon- DifferenceDifferenceDiffer- (m) 
ent ent (m/sec) (m/sec) ence 
Wind Wind (m) 
Differ- Differ- 
ence ence 


(m/sec) (m/sec) 


Persis- 

tence 4.19 5.45 TO OV 25 Sl.) 63.40 
None 4.82 4.35 7205 1, 29 106 (0s HAL Steve 
1 4.46 4.15 6.48 6.54 49.29 68.89 
2 4.32 G07, 6.24 6.78 48.57 O73 
3 a | 4.04 6613 67.0 48.18 o7 2709 
4 4.27 GyO2 oy ey JET8) 6.65 47.97 66.70 
5 4.30 4.03 oye Ey 6.63 NIP ES) 66.60 
6 4.34 4.04 9) s dbl 6.64 48.02 66.66 
7 4.40 4.04 6nZz5 6.64 48.15 CO. 75 
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Ertterineg up to four or five times improves the verifica-— 
tion values. However, beyond five filterings, the verifi- 
cation values begin to decline. The u wind component veri- 
fication is different, however. Area average modulus u 
Semponent differences remainhighner than persistence values. 
RMS u component differences with one pass of the filter are 
smaller than for a persistence forecast and continue to 
improve up to four or five passes. This can be explained by 
noting that wind gradients north Aa Sew of u component 
wind maximum are much larger than the gradients around v 
component wind maximum (see Figs. 14 and 16). Therefore, 

a slight displacement of the predicted maximum north or 
south of the actual position leads to larger errors than are 
produced in the v component by east-west phase errors. 

Even though filtering reduces the actual jet core speed, 

the diffusion of the maximum due to filtering results in a 
better verification value. 

Figs. 14 through 19 show the fields used for NVA prog- 
nosis verification, and the prognosis error fields for the 
v component, u component and height fields respectively. 
All fields, especially the v component and height field 
show that horizontal advection of most of the synoptic 
features was too rapid. Predicted positions and intensity 
of the major synoptic features are shown by the stars and 
values in parentheses on the verification charts. One can 
Pee that the analysed error values are predominantly a 


function of position, and not intensity or addition/deletion 
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Fig. 14. 300 mb NVA balanced v component wind verification 
analysis from 60E to 150W for 09/1200Z January 
1972. Star denotes prognosis wind maxima location 
With intensity in parenthesis. Units = m/sec. 








mipee to. Nondivergent input NVA prediction verification 
error of v component wind. Positive (+) value 
indicates too strong southerly wind. Units = 
m/sec. 
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Fig. 16. 300 mb NVA balanced u component wind verification 
analysis from 60E to 150W for 09/1200Z January 


1972. Star denotes prognosis wind maxima location 


with intensity in parenthesis. Units = m/sec. 
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Fig. 17. Nondivergent input NVA prediction verification 


error of u component wind. Positive (+) value 
indicates too strong westerly wind. Units = 
m/sec. 
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figs lo. S00 mb NVA balanced height verification field from 


60E to 150W for 09/1200Z January 1972. Trough 
Symbol denotes prognosis trough position, Units = 
decameters. 
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error of height field.) Units = m. 


a2 





of systems. A typical example is cir short wave trough 
southwest of Korea. The prediction model overpredicts the 
trough movement (See Fig. 18) with a resulting positive 
height error pattern west of the predicted trough position 
and a negative center east of the trough (see Fig. 19). 

The model moves both of the v component maxima associated 
with the trough too far south as shown in Fig. 14. This 
position error contributes about one-half of the value of 
the verification error analysis shown in Fig. 15. The re- 
Mainder of the error is due to an underestimate of the 
magnitude of the v component maxima, which seems to be true 
of the majority of the centers. The u component wind 
minima associated with the trough is forecast to move more 
rapidly eastward than the verifying position shown in Fig. 
16. This is depicted in Fig. 17 as a negative error center 
east of Korea and a positive center west of Korea. Unlike 
the v component winds, the prediction appears to retain u 


component magnitudes comparable with the verification. 


C. RESTORATION TECHNIQUES 

Results of the experiments with climatological input 
data had suggested that restoration of both wind fields and 
the height field at each time step also tended to generate 
excessive 2-d gravity waves originating from the boundaries. 
This can be explained as follows: As the wind and height 
fields in the region of outflow try to adjust to the fixed 


boundary values, the restoration coefficient allows only 
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partial adjustments. This produces a continual imbalance, 
which is manifested by gravity waves propagating from this 
region. The generation of these unwanted gravity waves 
can be reduced by either restoring the wind values and not 
the heights, or vice-versa. Slightly smaller amplitude 
gravity waves were produced, and better verification values 
were obtained, when u and v component winds were restored 
with no restoration of the height field. This result is 
shown in Table 5A for both the enol Staggered grid re- 
sults where Forecast value is defined as in Eq. 10 and the 
value in parenthesis is the ratio of RMS height error to a 


RMS height error from a persistence forecast. 


D. DIVERGENCE VALUES OF INPUT NVA FIELDS 

As mentioned earlier, a limitation of the divergence in 
the NVA input to the prediction model appeared to have de- 
Sirable effects on the elimination of unwanted two-grid- 
length gravity waves. This might be expected to have an 
effect on the verification values of the prediction, and is 
shown in Table 5B. The best verification was obtained using 
an NVA output made nondivergent input to the prediction 
(Method 2 in Table 3). Slightly worse results were obtained 
from NVA input data obtained from a nondivergent input into 
the NVA program (Method 1 in Table 3). Considerably less 
skill was achieved using initial data from the NVA program 
Mindi vergent input (Method 2a in Table 3). In fact, the 
full grid prediction model blew up before a 24-hour forecast 


could be obtained. One would expect a wind field which was 
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TABLE 5 


PROGNOSIS VERIFICATION AS A FUNCTION OF 


RESTORATION, 


Input 


A. RESTORATION 
Full Grid 


Restoration 
Restoration 


Staggered Grid 


Restoration 
Restoration 
Restoration 


NVA(1,7,25)made 
NVA(1,7,25)made 
NVA(1,7,25)made 


nondivergent 
nondivergent 
nondivergent 


B. NVA DIVERGENCE 


Full Grid 


NVA(1,7,25)with 
NVA(1,7,25) made 


Staggered Grid 


NVA(1,7,25)with 
ive (l,/7,25)made 


nondivergent 
nondivergent 


nondivergent 
nondivergent 


C. NVA WEIGHTING 


Full Grid 
Mmeecl.,1=1/2,25) 
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nondivergent to produce fewer gravity waves than a divergent 
wind field, and the resulting gravity waves appear to be one 
PmactOr responsible for degrading the verification results in 


the divergent NVA input wind fields. 


E. EFFECT OF NVA WEIGHTING 

A variation of the weighting factors in the NVA analysis 
could be expected to vary the mead = verification. A 
closer specification of the wind changes allowed in the NVA 
analysis (larger G@) should produce a better forecast wind 
field, as described in the diagnostic results section on 
NVA weighting. Table 5c shows this result for the stag- 
gered grid prediction. In the full grid, specifying a=25 
produced a computationally stable prediction for 24-hours, 
but using B = 7 or 1.5 did not adequately limit the changes 
in the wind components and the predictions became unstable 


before 24-hours. 


F. STAGGERED VS. NON-STAGGERED GRID 

Tables 5 and 6 point out that the staggered grid predic- 
tions in all cases produce more accurate forecasts than a 
full grid with comparable input. It is difficult to deter- 
mine if this is a result of the removal of the computational 
mode from the walls by use of the staggered grid or an even 
more Se€vere overprediction of the movement of systems in the 
full grid compared to the staggered grid. The excessive 
amplitude of the 2-d gravity waves evident in the full grid 


pradictions, as compared to the staggered grid, appears to 
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be the most significant influence. 


G. COMPARISON OF INITIALIZATION METHODS 

On both the staggered and the nonstaggered grids, the 
NVA fields made nondivergent input to the prediction had 
the best verification value. Considerably poorer were the 
balance restoration and balance no-flux, with the forecasts 
with restoration being slightly eee Table 6). 
Better verification of the NVA input might be expected 
Since the changes necessary in the diagnostic phase of 
prediction are smaller for winds and much smaller for 
height than in the nonlinear balance method. This results 
in individual synoptic features retaining their identity 


and as a result are forecast more accurately. 
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TABLE 6 


PROGNOSIS VERIFICATION AS A FUNCTION OF 
INITIALIZATION TYPE 


Input Restore Forecast 
Value 

Pall Grid 

No-Flux -20 + (.57) 

Restoration u & Vv 2 Suet (Olen) 

NVA(1,7,25)made nondivergent u & V 4.52 + (.84) 


staggered Grid 


No-Flux 3455 6.55) 
Restoration u & Vv 3.77 + (€.58) 
NVA(1,7,25)made nondivergent u & Vv 5.86 + (.95) 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


The staggered grid forecasts produced better verifica- 
tion statistics than the full grid. This is desirable since 
the barotropic full grid prediction requires approximately 
33 minutes of computer time for a 24-hour forecast, while 
the staggered grid requires only about 10 minutes. This is 
a considerable saving of time when poner des ie operational 
requirements. 

Restoration near the boundaries of the u and v component 
Winds and no restoration of the height field produced better 
verification values than restoring all fields or restoring 
only the height field. This results in a slight additional 
Peauction in computer time. 

A nondivergent NVA input to the prediction model pro- 
duced better verification statistics than either the balance 
restoration or balance no-flux technique. Nondivergent NVA 
input also produced better verification statistics than NVA 
fields which had no constraint on divergence. The smallest 
amount of adjusting of wind speed from the input fields, 
offset by rapidly increasing NVA convergence time, produced 
better verification values. To obtain these NVA fields re- 
quired running two programs with a combined run time of 
about thirty-five minutes. Again, this is excessive for an 
operational product. Three changes are possible to reduce 


the computation time. First, the NVA program could be 
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altered to place an additional constraint on the magnitude 
of the divergence of the resulting fields. This would 
eliminate the need for running the nondivergence progran. 
Also, further experiments might show that a specified 
amount of divergence results in a better initial field than 
no divergence. Also, time might be reduced by altering the 
values of the NVA coefficients to require less adherence to 
the input fields, and possibly larger tolerances on the 
horizontal acceleration. The most significant time saving 
could be realized by converting the NVA program to operate 
on a staggered grid. This would allow for a solution con- 
sidering only one-quarter of the present points. Since 
iterative processes are inherently more efficient on a 
smaller grid, the time savings should be well in excess of 
peractkor Of four. 

Boundary conditions on the Mercator belt still remain 
a problem. Research by Lieutenant Commander E. Harrison 
(personal communication) on appropriate boundary conditions 
for numerical models has suggested a different method of 
handling the prediction in regions of outflow. 

It is to be expected that a barotropic model, which can 
only advect synoptic features with the wind at one level, 
should produce errors in displacement; specifically a 300 mb 
level should overforecast movement. In the real atmosphere 
the vertical differences in both temperature and momentum 
advection can make the barotropic approximation a poor one. 


This would suggest that a baroclinic model would give 
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significantly better movement forecasts. The staggered grid 
now makes a three- or four-level model operationally feasi- 
ple. Initial experiments, on a s imiieed grid used by Harri- 
son (1969) and Steinbruck (1971), have shown that at least 

a three-level model is desirable. At present, the time step 
in the staggered prediction models is five minutes. This 
could be doubled if the advection of the height field in 

Eq. 5 was finite-differenced over twice the present space 
increment in the full grid, as is the case for momentum ad- 
vection, with a staggered grid in the time sense. The use 
of semi-implicit time differencing in the prediction models 
would also be a means for increasing the time step. One 
possible approach is a time-averaging form suggested by 
Brown (1971). 

By using at least three layers, it is possible to 
parameterize the effect of cumulus convection, which is so 
important in the tropics. A baroclinic model would also 
make experiments with tropical cyclones movement and 


development feasible. 
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